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Abstract 

The spatial structure of fluctuations in spatially inhomogeneous processes can be modeled in 
d ' terms of Gibbs random fields. A local low energy estimator (LLEE) is proposed for the interpo- 



lation (prediction) of such processes at points where observations are not available. The LLEE 



> 

! approximates the spatial dependence of the data and the unknown values at the estimation points 

• ^ ! by low-lying excitations of a suitable energy functional. It is shown that the LLEE is a linear, 

unbiased, non-exact estimator. In addition, an expression for the uncertainty (standard deviation) 

^ ; 

, of the estimate is derived. 
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I. INTRODUCTION 

Spatial random fields (SRF's) have applications in hydrology |9ljllj|. oil reservoir engi- 
neering g e_ental pollutant .nappmg and n.k a=_tJI„tWng exploration 

and reserves estimation |3], as well as environmental health studies SRF's model spatial 
correlations in variables such as mineral concentrations, dispersion of environmental pol- 
lutants, soil and rock permeability, and flow fields in oil reservoirs. Knowledge of spatial 
correlations enables (i) generating predictive iso-level contour maps (ii) estimating the uncer- 
tainty of predictions and (iii) developing simulations that partially reconstruct the process of 
interest. Geostatistics provides mathematical tools for these tasks. The classical approach 
is based on Gaussian SRF's (GSRF's) and various generalizations for non-Gaussian distri- 



butions 
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, For GSRF's the spatial structure is determined from the covariance matrix, 
which is estimated from the distribution of the data in space. 

An SRF state (realization) can be decomposed into a deterministic trend m^{s), a cor- 
related fluctuation X\{s), and an independent random noise term, e(s) so that X{s) = 
m^{s) + Xa(s) + e(s). The trend represents large-scale variations of the field, which can be 
obtained in principle by ensemble averaging, i.e. mx(s) = i?[X(s)]. In practice, the trend 
is often determined from a single available realization. The fluctuation term corresponds 
to 'fast variations' that reveal structure at small scales, which nonetheless exceed a cut-off 
A. The random noise represents non-resolved inherent variability due to resolution limits, 
purely random additive noise, or non-systematic measurement errors. It is typically assumed 
that the fluctuation is a second-order stationary SRF, or an intrinsic SRF with second-order 
stationary increments The observed SRF after detrending is a zero- mean fluctuation: 
X*{s)=Xx{s) + e{s). 

In statistical physics the probability density function (pdf) of a fluctuation field x(s) 
governed by an energy functional H[x{s)] is expressed as /x[x(s)] = exp {—H[x{s)]} , 
where Z is the partition function. Using this representation, the Gaussian joint pdf in 
classical geostatistics is expressed in terms of the functional: 

H[xis)] = ^jdsj ds'xis) [G^]-\s - s')x(s'). (1) 

In Eq. (^, [G'2;]^^(s — s') is the inverse of the covariance function Gx{s — s'), which determines 
the spatial disorder. While statistical physics plays an increasingly important role in under- 
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standing the behavior of complex geophysical systems its applications in geostatistical 
analysis have not yet been explored. 

Spartan Spatial Random Fields (SSK 
actions', in the spirit of Markov SRF's 



s) model spatial correlations in terms of 'inter- 

n 



14| . In 6| general properties and permissibility 
conditions were derived for the fluctuation-gradient-curvature (FGC) SSRF model, with the 
following energy functional: 

^^fgc[X,] = ^ I ds{[X,{s)f + r^,eNX,{s)f + e[V'X,{s)f]. (2) 

For this model, a moment-based method for parameter estimation was proposed and 
tested with simulated data; methods for SSRF non-constrained simulation were presented 



in |2|; systematic reduction of anisotropic disorder, based on the covariance tensor identity, 
was investigated in main parameters: the scale coefficient 

rjQ, the covariance-shape coefficient rji, and the correlation length ^. Bochner's theorem ^| 
for the covariance function requires rji > —2. A coarse-graining kernel is used to cut off 
the fluctuations at kc oc UQ, leading to band-limited covariance spectral density and 
differentiable field configurations (in the mean square sense) Q. 



II. OPERATOR NOTATION 

Let neR'^ denote the area of interest and A{Q) its boundary. Consider an SSRF defined 
over this area with parameters i]Q,rii,^, with a finite variance cr^. Let us assume that it is 
possible to normalize the SSRF to unit variance by simply dividing the states with the 
standard deviation. Next, it is possible to express the pseudo-energy functional in terms of 
an operator notation notation as follows: 

H[Xx] = {Xx{s)\n |Xa(s) ) + S{A) = [ dsX^is) H [Xa(s)] + S{A), (3) 

Jn 

where 7i is a 'pseudo-hamiltonian' operator and S{A) is a surface term. Assuming that the 
surface term is negligible, the eigenvalue equation becomes: 

H\M^-M) = EM^;^), (4) 

where ipE{s',h) is an eigenfunction, E is the corresponding energy and b a degeneracy 
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vector index, which may include both discrete and continuous components. Since the SSRF 
has been normahzed to unit variance, the eigenf unctions '?/'£;(s;b) can also be assumed 
normalized, i.e., dsip'^ls; b) = 1, and then -^[Xa] = E. 

If Eq. (0]) admits solutions for non-zero E, one can construct eigenfunctions that cor- 
respond to positive excitation energies E. The realization probability that corresponds to 
low-lying excitations is high. Hence, the main idea is to consider the observed state or the 
union of the observations and the predictions as being locally represented by an excited 
state. This approach can be used for both parameter estimation and prediction (spatial 
estimation) 

A. Eigenfunctions for FGC case 

For the FGC functional of Eq. integrating the square-gradient term by parts leads 
to the following equation: 

/ ds[VM^-M^ = - I rfs^i5(s;b)V2^ij(s;b)+ / rfa ■ V^ij(s; b) ^^^(s; b). (5) 
Jn Jn JA(n) 

In Eq. (0) /4(f^) da denotes the surface integral on the boundary of the area of interest. 

Secondly, using Green's theorem on the square- curvature term one obtains 

/ ds [V^ijE{s;h)Y = [ tis^s(s;b)VVi,(s;b)+ / da ■ WM^;h)V''M^-M 
Jn Jn J A{n) 

- f da-V[VVi?(s;b)]^E(s;b). (6) 

JA{n) 

Hence, in the operator notation the FGC functional is expressed as follows: 

^%c = ^[l-r/ie^V^ + e^V^], (7) 
and the surface term is given by: 



s{n) 



?7i^^ / da-Vi)E{s;h)'ipE{s;h) 
JA(n) 



+ / rfa ■ VV^i,(s;b)VVs(s;b) 

JA{n) 

-el da-V[V^^E{s-M]^E{^M 
J Am 



(8) 



If the units are chosen so that 2t]q^'^ = 1 and the surface term is ignored, the eigenvalue 
equation is given by the following partial differential equation (pde): 

7/>s(s; b) - r/i e V ^e{s; b) + ^ ^e{s; h) = E ^^(s; b). (9) 
The eigenfunctions ipEi^', b) of Eq. are given by the following four plane waves: 

V'i5(s;b) = e*^^", kj = kjO, (10) 

where represents the unit direction vector, and kj the magnitudes of the characteristic 
wave-vectors that are given by the roots of the fourth- order characteristic polynomial: 

nfg,(A;0 = il-E) -7^^ e + ^ = 0. (11) 

Thus, the characteristic wavevectors are given by the following expressions: 
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k,{vu^,E) = -^^\jr^, + ^lr^l-A{l-E) (12) 
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k,{7i,,i,E) = --^^\lr^, + ^/r^l-A{l-E) (13) 
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k,{v,,^,E) = -^^\lr^,-^lr,l-A{l-E) (14) 
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fc4(r/i,e,i?) = -^^\jvi-\jril-Kl-E). (15) 

Note that only the magnitude of the wave- vectors is determined from the pde Q. This is 
due to the fact that isotropic spatial dependence was assumed in the SSRF model. 

(a) If ?7i > A 1 - r7^/4 < E < 1 all the roots are real, (b) If ?7i > A E > 1, then ki, k2 
are real, while /ca, k^^ are purely imaginary, (c) If r/i > A 1 — ?7^/4 > i?, then all the roots 
are complex, (d) If r^i < A 1 — 7/^/4 < E < 1, then all the roots are imaginary, (e) If 
Tji < OAE > 1, then ki, k2 are real, while k^, k^ are imaginary, (f) If rji < OAO < E < l—rjl/i 
all the roots are complex. 

In general, an excited state formed by the linear superposition of degenerate eigenstates 
of energy E is given by the expression: 

ZEis; Cb) = Y,u{k,- II A;,- II) / decj{e) exp (^k, ^ ■ s) , (16) 
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where Cj{6) is a direction-dependent (possibly complex- valued) function, \\kj\\ is the modulus 
of the characteristic wavevector, and u{.) is the unit step function, used to guarantee that the 
fluctuations in the excited state do not exceed the cutoff 'frequency'. For the estimation of 
real- valued processes, the coefficients Cj{0) are constrained to give real values for the excited 
state Ze{s] Cb). If Cj{6) = Cj, an isotropic excited state is obtained, which can be expressed 
as Ze{s; Ci, . . . ,04) = Yl'^j=iCj u {k^ - ||/Cj||) ^^(s; j), where i^E{s;j) = f dOexp {kj 6 ■ s^. 

B. Eigenstates in d = 1 

We examine in more detail the real- valued eigenstates that are trigonometric or hyperbolic 
functions in the one-dimensional domain [0, L] G M. 

1. Exponential Eigenstates 

For characteristic wave-vectors k that are real numbers, the normalized eigenf unctions 
and the corresponding energies of Eq. are given by 

B = + (18) 

However, if the exponential function is inserted in Eq. 0, the resulting energy is given by 

H[Xis)] = l + ViikO' + ikO'- (19) 

The difference between the energy given by Eq. (fTTj) and the correct energy, given by Eq. (fTUj) 
is due to the fact that the boundary term can not be ignored for the localized exponential 
excitation. 

2. Trigonometric Eigenstates 

If k is an imaginary number, the eigenfunctions are trigonometric functions. A normalized 
cosine eigenfunction and the corresponding energy are given by: 
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For large domains, kL » 1, Eq. (PT|) is practically equivalent to Eq. (fT^. As expected, in 
the case of an extended eigenstate (as the cosine) the boundary term can be ignored. 



III. SPATIAL ESTIMATION WITH SSRF'S 



Assume 5*111 = (si, • • -Sjv) is a set of data points with the respective vector of measure- 
ments denoted by X* = (Xj*, . . . ,X^); let Sq ^ be the estimation point and Xa(so) the 
estimate (spatial prediction). The local neighborhood of Sq is the set So = -B(so; Tc) of all 
the data points Sj,j = 1,...,M inside a 'sphere' of radius equal to one correlation range 
:'rom Sq. In geostatistics, X(so) is determined by optimal linear filters (kriging estimators) 
^, Q , which form the estimate as a superposition of the data values inside the local neigh- 
borhood, and there is no explicit resolution scale. The coefficients of the superposition are 
selected to make the estimate unbiased and to minimize the mean square error. Kriging 
is an exact interpolator, meaning that for any Sj G Sni,X{si) = X*(sj). Exactitude is not 
always desirable, since it ignores measurement errors and leads to excessive smoothing of 
the fluctuations. Hence, different estimation methods are useful. The SSRF models can be 
used in kriging algorithms to provide new, differentiable covariance functions. In addition, 
within the SSRF framework it is possible to deflne a new type of estimator. 



A. Low Local Energy Estimators 

The central idea is that a 'good' estimate should correspond to a state with signifl- 
cant probability of realization. If the energy functional is non- negative, as in Eq. (jSj), the 
highest probability is associated with the uniform state Xa(s) = 0, which is not phys- 
ically interesting. Other states with high probability correspond to low-energy excita- 
tions. Let us superimpose the degenerate eigenstates with energy E to form a mixed state 
Ze{s] c) = YliLi V'£;(s; &«); c = (ci, . . . , cd) is a D-dimensional vector of linear coefficients 
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that correspond to the degeneracy indices. In principle D can be infinite since the direc- 
tional dependence given by Eq. (fTBj) is continuous. However, in practice it may be simplest 
to restrict the search to one 'optimal' direction. The energy H[Ze{s; Cb)] of the mixed state 
is not necessarily equal to E. In fact, for orthonormal eigenstates H[Ze{s; c)] = fiE, where 
fj, = J2f=i ^i- This refiects the fact that the 'energy level' of the observed process is set by the 
measurements (i.e., the coefficients q). Since the scale coefficient r^o is inversely proportional 
to the magnitude of the fluctuations, it follows that oc i]q. It should also be noted that 
if two mixed states (ci,i?i) and (02,-^2) are energetically equivalent, i.e., fiiEi = /i2-E'2, 
they are not in general linearly related, since according to Eqs. (fTn|l . (|T^ - (fTH|l and (|TH|l . the 
dependence of Ze{s; c) on E is nonlinear. 

We propose that the observations for G -B(so; r^) be expressed as X*{sj) = ZE^Sj] Cq) + 
e{sj), where Ze{s]Cq) is a 'local' excitation and £{sj) is the local excitation residual. Local 
dependence stems from the fact that the coefficients Cq depend on Sq, in contrast with the 
solution of Eq. (fT^ . in which the coefficient vector is global. The LLEE estimator is then 
given by Xa(so) = ^_e(so; Cq). Since Ze{p] Cq) is an estimate of the underlying process Xa(s), 
the excitation residual s{sj) is not in general the same as the noise e(s). The coefficients Cq, 
follow from minimizing the mean square excitation residual inside -B(so; r^) , i.e.. 



The above is a typical problem of multiple linear regression, where the regressors are the 
functions '>Pe{^u bj)- If we define the M x D matrix i'E,ij = i^Ei^i'jbj), the solutions for co,i 
and the LLEE are given by: 



M 




c 



(22) 



M 



'E,jk, 



I = 



1,...,M; k = l,...,D, 



(23) 



D 



M 




(24) 



k=l 



1=1 



Xa(so) = Wo • X*, 



(25) 
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where wq is a weight vector given by: 

D D 

Wo,i = ^ ipEflk ^ [a]k} ^E,ij, i = l,...,M. (26) 

k=i j=i 

The uncertainty of the LLEE estimate is determined from the ensemble variance of the 



local excitation residual (t'^{sq) = E 



X*(so)-X(so) 



I.e.: 



,0i 



(27) 



MM M 

O-g^(So) = f^x* + X] X] *-^^*'*^' ~ ^ 2Z ^0'^ '-^^* 

i=l j=l j=l 

where G^.^ij = E[X*X*] is the covariance matrix at the observation points, Gx*fli — 
E [X* Xq], is the covariance vector of the fluctuations between Sq and the estimation point, 
and cr^, = E [X* X*] is the variance of the observed process. 



B. Properties of the LLEE 

It follows from Eqs. and ()26|) that the LLEE is linear in the fluctuations. Hence, 
the estimates are unbiased and follow the Gaussian law (if the observations are normally 
distributed). Kriging methods are based on minimization of the (ensemble) mean square 
error, which is a global optimality criterion. In contrast, the LLEE criterion is local (i.e., 
minimum of the average squared excitation residual in the neighbourhood of the estimation 
point). Another difference with kriging is that low local energy estimates do not match 
exactly the measurements at observation points. The property of non-exactitude is main- 
tained even when the noise can be ignored. Finally, unlike kriging predictions, the LLEE 
provides multiple estimates, since different energy levels lead to different excited states. In 
this respect the LLEE is similar to a simulation method. However, simulations involve the 
generation of random numbers, in contrast with the LLEE method. It should also be noted 
that the energy of local excitations is not necessarily the energy of the estimated state, be- 
cause the locality of the coefficient vector Cq means that the operators V and contribute 
to the overall energy when they act on the coefficients of the mixed state in Eq. Q. 



IV. CONCLUSIONS 



A spatial estimation method for applications in the geosciences is presented. The method 
is based on the use of 'pseudo-energy' functional, motivated by exphcit constraints or 
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heuristic physical arguments, to capture the spatial heterogeneity of the observed process. 
Estimates of the process at unmeasured points (predictions) are based on local interpolating 
functions that represent low-energy excitations of the pseudo-energy. Multiple estimates of 
the process can be generated by considering local interpolating functions that correspond to 
different excitation energies. 
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